#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <fstream>
#include <math.h>
#include <cstdlib>
#include <unistd.h>
#include <glob.h>
#include <vector>

using namespace std;
#include "asteroid.cpp"
#include "asteroidpelem.cpp"
#include "asteroidwithsize.cpp"
#include "asteroidwithtax.cpp"
#include "asteroidwithpole.cpp"
#include "masiero.cpp"
#include "nugent.cpp"
#include "simps.cpp"
#include "akari.cpp"
#include "alilagoa.cpp"
#include "carvano.cpp"
#include "demeocarry.cpp"
#include "cibuki.cpp"

long countFileLines(string fileName)
{
  ifstream myfile (&fileName.at(0));
  long i=0;
  string line;
  
  if (myfile.is_open())
  {    
    while ( getline (myfile,line))
      if (line.length()>1)
	i++;
    myfile.close();
  }
  else cout << "Unable to open file";
  return i;
}

void countMPCorb(long *nMPC, string fileName)
{
  ifstream myfile (&fileName.at(0));
  long i=0;
  string line;

  *nMPC=0;
  
  if (myfile.is_open())
  {
    for (int j=0; j<41; j++)  // MPC file has 
      getline (myfile,line) ; // skip header
    
    while ( getline (myfile,line))
      if (line.length()>1)
	  i++;
    *nMPC=i;
    myfile.close();
  }
  else cout << "Unable to open file"; 
}

template<class T>
void readMPCorb(T *myAst, long *indLastMPC, long *indLastNumbMPC, string fileName)
{
  ifstream myfile (&fileName.at(0));
  long i=0;
  string line;

  *indLastNumbMPC=0;
  *indLastMPC=0;
  
  if (myfile.is_open())
  {
    for (int j=0; j<41; j++)  // MPC file has 
      getline (myfile,line) ; // skip header
    
    while ( getline (myfile,line))
      if (line.length()>1)
	{
	  myAst[i].MPCstring2elem(line);
	  if ( (myAst[i].Number<=0) && (*indLastNumbMPC==0))
	    *indLastNumbMPC=i-1; 
	  //	  printf("%ld \t ",i); myAst[i].print();
	  i++;
	}
    *indLastMPC=i-1;
    myfile.close();
  }
  else cout << "Unable to open file"; 

  /*
  cout << "Last Numberd Index " << *indLastNumbMPC << "\n";
  myAst[*indLastNumbMPC].print();
  cout << " Last index " << *indLastMPC << "\n" ;
  myAst[*indLastMPC].print();
  */
}

template<class T>
void readMPCorb(T *myAst, long *indLastMPC, long *indLastNumbMPC)
{
  char fileName[]="MPCORB.DAT";
  readMPCorb(myAst, indLastMPC, indLastNumbMPC, fileName);
}

// attempt to find ast within the arr UN-NUMBERED
// the array must be sorted
template<class T>
long findUnnumberedAst(T *arr, T ast, long lastInd, long firstInd=0)
{
  long ini=firstInd;
  long end=lastInd;
  long pivot = (end+ini)/2;

  //  cout << "Init " << ini << " ";  arr[ini].print();
  //  cout << "Last " << end << " ";  arr[end].print();
  //  cout << "Ref  ";  ast.print();

  long iter=0;
  
  while ((pivot-ini)>1)
    {
      if (ast==arr[pivot]) return pivot;
      if ((ast<arr[pivot]))
	  end=pivot;
      else
	  ini=pivot;
      pivot = (end+ini)/2;
      // cout << "iter= " << iter++ << " init= " << ini << " end= " << end << " pivot= " << pivot << " "; arr[pivot].print();
    }
  
  for(long j=ini; j<=end; j++)
    if (ast==arr[j]) return j;
  
  return -1;
}

template<class T>
void readInAllSYN(T myAst[], long indLastMPC, long indLastNumbMPC, string* ids, long nIds)
{
  string line;
  ifstream myfile ("all.syn");
  long ind;
  AsteroidWithPole ast;

  if (myfile.is_open())
  {
    for (int j=0; j<2; j++)
      getline (myfile,line) ; // skip header
    
    while ( getline (myfile,line))
      if (line.length()>1)
	{
	  ast.SYNstring2propElem(line);
	  if (ast.Number>0)
	    {
	      myAst[ast.Number-1].ap=ast.ap;
	      myAst[ast.Number-1].ep=ast.ep;
	      myAst[ast.Number-1].sinip=ast.sinip;
	    }
	  else
	    {
	      // can you match is right there ?
	      if ((ind = findUnnumberedAst(myAst,ast,indLastMPC,indLastNumbMPC+1))>=0) // YES ! 
		{
		  myAst[ind].ap=ast.ap;
		  myAst[ind].ep=ast.ep;
		  myAst[ind].sinip=ast.sinip;
		  //		  cout << ">>match 1  " << ast.pack << endl;
		}
	      else
		{ // the designation could not be matched.
		  // try to see if the asteroid has been numbered LATER
		  ast.findId(ids,nIds);
		  if (ast.Number>0)
		    {
		      myAst[ast.Number-1].ap=ast.ap;	  
		      myAst[ast.Number-1].ep=ast.ep;	  
		      myAst[ast.Number-1].sinip=ast.sinip;
		      //		      cout << ">>match 2  " << ast.pack << endl;
		    }
		  else
		    {
		      if ((ind = findUnnumberedAst(myAst,ast,indLastMPC,indLastNumbMPC+1))>=0) // YES ! 
			{
			  myAst[ind].ap=ast.ap;	     
			  myAst[ind].ep=ast.ep;	     
			  myAst[ind].sinip=ast.sinip;
			  //			  cout << ">>match 23  " << ast.pack << endl;
			}
		      else
			{
			  cout << ">>match FAILED  " << ast.pack << "--" << line << endl;
			}
		    }
		}
	    }
	}
    myfile.close();
  }
  else cout << "Unable to open file"; 
}

template<class T>
void readInUfitANA(T myAst[], long indLastMPC, long indLastNumbMPC, string* ids, long nIds)
{
  string line;
  ifstream myfile ("ufitobs.pro");
  long ind;
  AsteroidWithPole ast;

  if (myfile.is_open())
  {
    for (int j=0; j<3; j++)
      getline (myfile,line) ; // skip header
    
    while ( getline (myfile,line))
      if (line.length()>1)
	{
	  ast.uFITstring2propElem(line);
	  if (ast.Number>0)
	    {
	      myAst[ast.Number-1].ap=ast.ap;
	      myAst[ast.Number-1].ep=ast.ep;
	      myAst[ast.Number-1].sinip=ast.sinip;
	    }
	  else
	    {
	      // can you match is right there ?
	      if ((ind = findUnnumberedAst(myAst,ast,indLastMPC,indLastNumbMPC+1))>=0) // YES ! 
		{
		  myAst[ind].ap=ast.ap;
		  myAst[ind].ep=ast.ep;
		  myAst[ind].sinip=ast.sinip;
		  //		  cout << ">>match 1  " << ast.pack << endl;
		}
	      else
		{ // the designation could not be matched.
		  // try to see if the asteroid has been numbered LATER
		  ast.findId(ids,nIds);
		  if (ast.Number>0)
		    {
		      myAst[ast.Number-1].ap=ast.ap;	  
		      myAst[ast.Number-1].ep=ast.ep;	  
		      myAst[ast.Number-1].sinip=ast.sinip;
		      //		      cout << ">>match 2  " << ast.pack << endl;
		    }
		  else
		    {
		      if ((ind = findUnnumberedAst(myAst,ast,indLastMPC,indLastNumbMPC+1))>=0) // YES ! 
			{
			  myAst[ind].ap=ast.ap;	     
			  myAst[ind].ep=ast.ep;	     
			  myAst[ind].sinip=ast.sinip;
			  //			  cout << ">>match 23  " << ast.pack << endl;
			}
		      else
			{
			  cout << ">>match FAILED  " << ast.pack << "--" << line << endl;
			}
		    }
		}
	    }
	}
    myfile.close();
  }
  else cout << "Unable to open file"; 
}

template<class T>
void readInAllNumANA(T myAst[], long indLastMPC, long indLastNumbMPC, string* ids, long nIds)
{
  string line;
  ifstream myfile ("allnum.pro");
  long ind;
  AsteroidWithPole ast;

  if (myfile.is_open())
  {
    for (int j=0; j<3; j++)
      getline (myfile,line) ; // skip header
    
    while ( getline (myfile,line))
      if (line.length()>1)
	{
	  ast.AllNUMstring2propElem(line);
	  if (ast.Number>0)
	    {
	      myAst[ast.Number-1].ap=ast.ap;
	      myAst[ast.Number-1].ep=ast.ep;
	      myAst[ast.Number-1].sinip=ast.sinip;
	    }
	  else
	    {
	      // can you match is right there ?
	      if ((ind = findUnnumberedAst(myAst,ast,indLastMPC,indLastNumbMPC+1))>=0) // YES ! 
		{
		  myAst[ind].ap=ast.ap;
		  myAst[ind].ep=ast.ep;
		  myAst[ind].sinip=ast.sinip;
		  //		  cout << ">>match 1  " << ast.pack << endl;
		}
	      else
		{ // the designation could not be matched.
		  // try to see if the asteroid has been numbered LATER
		  ast.findId(ids,nIds);
		  if (ast.Number>0)
		    {
		      myAst[ast.Number-1].ap=ast.ap;	  
		      myAst[ast.Number-1].ep=ast.ep;	  
		      myAst[ast.Number-1].sinip=ast.sinip;
		      //		      cout << ">>match 2  " << ast.pack << endl;
		    }
		  else
		    {
		      if ((ind = findUnnumberedAst(myAst,ast,indLastMPC,indLastNumbMPC+1))>=0) // YES ! 
			{
			  myAst[ind].ap=ast.ap;	     
			  myAst[ind].ep=ast.ep;	     
			  myAst[ind].sinip=ast.sinip;
			  //			  cout << ">>match 23  " << ast.pack << endl;
			}
		      else
			{
			  cout << ">>match FAILED  " << ast.pack << "--" << line << endl;
			}
		    }
		}
	    }
	  
	}
    myfile.close();
  }
  else cout << "Unable to open file"; 
}

template<class T>
void tagAsteroidsWithC51(T myAst[], long indLastMPC, long indLastNumbMPC, string* ids, long nIds)
{
  ifstream myfile ("list.NumObsPacked");
  long ind;
  string line;
  char tax[5];
  AsteroidWithPole ast;

  // tag the numbered asteroids ==============================
  // reading the list NumObsPacked ============== with C51 ===
  if (myfile.is_open())
  {
    while ( getline (myfile,line))
      {
	if (line.length()>1)
	  {
	    ast.reset();
	    ast.line2pack(line);
	    ast.MPC2DesignationCompact();
	    if ((ind=ast.match(myAst,indLastMPC,indLastNumbMPC+1, ids, nIds))>=0) {
	      myAst[ind].hasWiseObs=1;
	    }
	    else {
	      cout << ">>match FAILED  " << ind << " " << ast.Designation << " -- " << line << endl;
	    }
	  }
      }
    myfile.close();
  }
  cout << "Numbered Tagged with C51" << endl;
  
  // tag the unnumbered asteroids ==========================
  ifstream myfileUnNumb ("list.UnnObsPacked");
  if (myfileUnNumb.is_open())
  {
    while ( getline (myfileUnNumb,line))
      {
	if (line.length()>1)
	  {
	    ast.reset();
	    ast.line2pack(line);
	    ast.MPC2DesignationCompact();
	    if ((ind=ast.match(myAst,indLastMPC,indLastNumbMPC+1, ids, nIds))>=0) {
	      myAst[ind].hasWiseObs=1;
	    }
	    else {
	      cout << ">>match FAILED  " << ast.Designation << " -- " << line << endl;
	    }
	  }
      }
    myfileUnNumb.close();
  }
  else cout << " Error opening UnnObsPacked file" << endl; cout.flush();

  cout << "UnNumbered Tagged with C51" << endl;  
}

//scanning asteroid 2009SR359 K09SZ9R
//>>match FAILED  2009SR359--K09SZ9R
//scanning asteroid 2009SY359 K09SZ9Y


string* getIds(long *N)
{

  long NnumIds=countFileLines("numids.txt");
  long Nids=countFileLines("ids.txt");
  long Ndbl=countFileLines("dbl.txt");
  string* ids = new string[NnumIds+Nids+Ndbl];
  string line;

  long i=0;
  *N = NnumIds+Nids+Ndbl;
  
  ifstream myfile ("numids.txt");
  if (myfile.is_open())
    {    
      while ( getline (myfile,line))
	if (line.length()>1)
	  ids[i++]=line;
      myfile.close();
    }
  else cout << "Unable to open file";

  ifstream myfile2 ("ids.txt");
  if (myfile2.is_open())
    {    
      while ( getline (myfile2,line))
	if (line.length()>1)
	  ids[i++]=line;
      myfile2.close();
    }
  else cout << "Unable to open file";

  ifstream myfile3 ("dbl.txt");
  if (myfile3.is_open())
    {    
      while ( getline (myfile3,line))
	if (line.length()>1)
	  ids[i++]=line;
      myfile3.close();
    }
  else cout << "Unable to open file";

  return ids;
}

template <class T >
int ast2bin(T *myAst, long N)
{
  ofstream fout("ast.bin");
  fout.write((char*)&N,sizeof(N));
  for (long i=0; i<N; i++)
    fout.write((char*)&myAst[i],sizeof(myAst[0]));
  fout.close();
}

AsteroidWithPole* bin2ast(long *N, long *indLastMPC, long *indLastNumbMPC)
{
  ifstream fin("ast.bin");
  fin.read((char*)N,sizeof(*N));

  AsteroidWithPole* myAst = new AsteroidWithPole[*N];

  indLastNumbMPC[0]=0;
  indLastMPC[0]=0;
  
  for (long i=0; i<N[0]; i++)
    {
      fin.read((char*)&myAst[i],sizeof(myAst[0]));
      if (myAst[i].Number>0) indLastNumbMPC[0]++;
      indLastMPC[0]++;
    }
  fin.close();

  return myAst;
}


template<class T>
void readInFamiliesDavidN(T myAst[], long indLastMPC, long indLastNumbMPC, string* ids, long nIds)
{
  long astNumber;
  string line;
  char tax[5];
  AsteroidWithPole ast;
  glob_t glob_result;
  ifstream familyFile;
  //  long fileLineCount;
  long familyNumber=0;
  double Cj;
  long ind;
  float H, ap, ep, sinip;

  glob("familiesDavidN/*.tab",GLOB_TILDE,NULL,&glob_result);

  for(unsigned int i=0; i<glob_result.gl_pathc; ++i)
    {
      cout << "Including family " << glob_result.gl_pathv[i] << "\t\t";
      familyFile.open(glob_result.gl_pathv[i]);

      if (familyFile.is_open())
	{
	  //	  fileLineCount=0;
	  while ( getline (familyFile,line))
	    {
	      if (line.length()>1)
		{
		  //		  astNumber = atol(&line.at(0));
		  sscanf(&line.at(0),"%ld %f %f %f %f %lf %d %ld",
			 &astNumber, &ap, &ep, &sinip, &H, &Cj, &ind, &familyNumber);
		  if ((astNumber>0) && (familyNumber>0))
		    {
		      myAst[astNumber-1].familyDavidN=familyNumber;
		      myAst[astNumber-1].Cj = Cj;
		    }
		  /*
		  if ((astNumber>0) && (fileLineCount==0))
		    {
		      familyNumber=astNumber;
		      cout << " family number is " << familyNumber << endl;
		    }
		  */
		  if (astNumber>0)
		    myAst[astNumber-1].familyDavidN=familyNumber;
		  else
		    cout << "big error occcurred " << endl;
		  //		  fileLineCount++;
		}
	    }
	  familyFile.close();
	} else {
	cout << "Opening family file failed " << endl;
      }
    
    }

  return;
}

template<class T>
void readInRotationPeriods(T myAst[], long indLastMPC, long indLastNumbMPC, string* ids, long nIds)
{
  ifstream myfile ("LC_SUM_PUB.TXT");
  long ind;
  AsteroidWithPole ast;
  long Number;
  char dummyStr[20];
  string line;
 
  if (myfile.is_open())
  {
    for (int j=0; j<5; j++)
      getline (myfile,line) ; // skip header
    
    while ( getline (myfile,line))
      if (line.length()>1)
	{
	  ast.LC_sumString2Asteroid(line);
	  if (ast.Number>0)
	    {
	      myAst[ast.Number-1].rotationPeriod=ast.rotationPeriod;
	    }
	  else
	    {
	      // can you match is right there ?
	      if ((ind = findUnnumberedAst(myAst,ast,indLastMPC,indLastNumbMPC+1))>=0) // YES ! 
		{
		  myAst[ind].rotationPeriod=ast.rotationPeriod;
		}
	      else
		{ // the designation could not be matched.
		  // try to see if the asteroid has been numbered LATER
		  ast.findId(ids,nIds);
		  if (ast.Number>0)
		    {
		      myAst[ast.Number-1].rotationPeriod=ast.rotationPeriod;	  
		      //		      cout << ">>match 2  " << ast.pack << endl;
		    }
		  else
		    {
		      if ((ind = findUnnumberedAst(myAst,ast,indLastMPC,indLastNumbMPC+1))>=0) // YES ! 
			{
			  myAst[ind].rotationPeriod=ast.rotationPeriod;	     
			  //			  cout << ">>match 23  " << ast.pack << endl;
			}
		      else
			{
			  cout << ">>match FAILED  " << ast.pack << "--" << line << endl;
			}
		    }
		}
	    }

	}
  }
}

template<class T>
void readInEarnDLR(T myAst[], long indLastMPC, long indLastNumbMPC, string* ids, long nIds)
{
  ifstream myfile ("table1_new.html");
  long ind;
  string line;
  char tax[5];
  AsteroidWithPole ast;
  int isBin=0;
  
  if (myfile.is_open())
  {
    // skip header
    for (int i=0; i<14; i++)
      getline (myfile,line);
    
    while ( getline (myfile,line))
      if ((line.length()>1)&&(line.at(1)=='<'))
	{
	  if (line.at(0)=='B')
	    if (isBin==1)
	      {isBin=0; continue; }
	    else {isBin=1;}
	    
	  // convert linefile to asteroid
	  ast.EARNString2Asteroid(line);
	  //	  cout << ast.pack << " "<< ast.Designation << " " << ast.Number << endl;
      	  // if the asteroid is Numbered in the file assing it
	  if (ast.Number>0)
	    {
	      myAst[ast.Number-1].D=ast.D;
	      myAst[ast.Number-1].sigmaD=ast.sigmaD;
	      myAst[ast.Number-1].pV=ast.pV;
	      myAst[ast.Number-1].sigmapV=ast.sigmapV;
	      //	      myAst[ast.Number-1].refs[3]+=MASIERO11;
	    }
	  else
	    {
	      // can you match is right there ?
	      if ((ind = findUnnumberedAst(myAst,ast,indLastMPC,indLastNumbMPC+1))>=0) // YES ! 
		{
		  myAst[ind].D=ast.D;
		  myAst[ind].sigmaD=ast.sigmaD;
		  myAst[ind].pV=ast.pV;
		  myAst[ind].sigmapV=ast.sigmapV;
		  //		  myAst[ind].refs[3]+=MASIERO11;
		  //cout << ">>match 1  " << ast.pack << "  "<< ind << endl;
		}
	      else
		{ // the designation could not be matched.
		  // try to see if the asteroid has been numbered LATER
		  ast.findId(ids,nIds);
		  if (ast.Number>0)
		    {
		      myAst[ast.Number-1].D=ast.D;
		      myAst[ast.Number-1].sigmaD=ast.sigmaD;
		      myAst[ast.Number-1].pV=ast.pV;
		      myAst[ast.Number-1].sigmapV=ast.sigmapV;
		      //		      myAst[ast.Number-1].refs[3]+=MASIERO11;
		      //cout << ">>match 2  " << ast.pack << endl;
		    }
		  else
		    {
		      if ((ind = findUnnumberedAst(myAst,ast,indLastMPC,indLastNumbMPC+1))>=0) // YES ! 
			{
			  myAst[ind].D=ast.D;
			  myAst[ind].sigmaD=ast.sigmaD;
			  myAst[ind].pV=ast.pV;
			  myAst[ind].sigmapV=ast.sigmapV;
			  //			  myAst[ind].refs[3]+=MASIERO11;
			  //cout << ">>match 23  " << ast.pack << endl;
			}
		      else
			{
			  cout << ">>match FAILED  " << ast.pack << "--" << line << endl;
			}
		    }
		}
	    }
	}
    myfile.close();
  } else cout << "file result_MarcoTargets-2.txt not found " <<endl;
}


template<class T>
void readInNEASourceRegion(T myAst[], long indLastMPC, long indLastNumbMPC, string* ids, long nIds)
{
  ifstream myfile ("nea.source.aei.dat.idl");
  long ind;
  string line;
  char tax[5];
  AsteroidWithPole ast;
  int isBin=0;
  
  if (myfile.is_open())
  {    
    while ( getline (myfile,line))
      if ((line.length()>1))
	{	    
	  // convert linefile to asteroid
	  ast.SourceRegionString2Asteroid(line);
	  //	  cout << ast.pack << " "<< ast.Designation << " " << ast.Number << endl;
      	  // if the asteroid is Numbered in the file assing it
	  if (ast.Number>0)
	    {
	      memcpy(myAst[ast.Number-1].neaSourceRegion, ast.neaSourceRegion, 10*sizeof(ast.neaSourceRegion[0]));
	      cout << ">>match 0  " << ast.pack << "  "<< ind << endl;
	      //	      myAst[ast.Number-1].refs[3]+=MASIERO11;
	    }
	  else
	    {
	      // can you match is right there ?
	      if ((ind = findUnnumberedAst(myAst,ast,indLastMPC,indLastNumbMPC+1))>=0) // YES ! 
		{
		  memcpy(myAst[ind].neaSourceRegion, ast.neaSourceRegion, 10*sizeof(ast.neaSourceRegion[0]));
		  //		  myAst[ind].refs[3]+=MASIERO11;
		  cout << ">>match 1  " << ast.pack << "  "<< ind << endl;
		}
	      else
		{ // the designation could not be matched.
		  // try to see if the asteroid has been numbered LATER
		  ast.findId(ids,nIds);
		  if (ast.Number>0)
		    {
		      memcpy(myAst[ast.Number-1].neaSourceRegion, ast.neaSourceRegion, 10*sizeof(ast.neaSourceRegion[0]));
		      //		      myAst[ast.Number-1].refs[3]+=MASIERO11;
		      cout << ">>match 2  " << ast.pack << endl;
		    }
		  else
		    {
		      if ((ind = findUnnumberedAst(myAst,ast,indLastMPC,indLastNumbMPC+1))>=0) // YES ! 
			{
			  memcpy(myAst[ind].neaSourceRegion, ast.neaSourceRegion, 10*sizeof(ast.neaSourceRegion[0]));
			  //			  myAst[ind].refs[3]+=MASIERO11;
			  cout << ">>match 23  " << ast.pack << endl;
			}
		      else
			{
			  cout << ">>match FAILED  " << ast.pack << "--" << line << endl;
			}
		    }
		}
	    }
	}
    myfile.close();
  } else cout << "file nea.source.aei.dat.idl not found " <<endl;
}

template<class T>
void readInNEASourceRegionGranvik(T myAst[], long indLastMPC, long indLastNumbMPC, string* ids, long nIds)
{
  ifstream myfile ("newprob.txt");
  long ind;
  string line;
  char tax[5];
  AsteroidWithPole ast;
  int isBin=0;
  
  if (myfile.is_open())
  {    
    while ( getline (myfile,line))
      if ((line.length()>1))
	{	    
	  // convert linefile to asteroid
	  ast.SourceRegionGranvikString2Asteroid(line);
	  //	  cout << ast.pack << " "<< ast.Designation << " " << ast.Number << endl;
      	  // if the asteroid is Numbered in the file assing it
	  if (ast.Number>0)
	    {
	      memcpy(myAst[ast.Number-1].neaSourceRegion, ast.neaSourceRegion, 10*sizeof(ast.neaSourceRegion[0]));
	      cout << ">>match 0  " << ast.pack << "  "<< ind << endl;
	      //	      myAst[ast.Number-1].refs[3]+=MASIERO11;
	    }
	  else
	    {
	      // can you match is right there ?
	      if ((ind = findUnnumberedAst(myAst,ast,indLastMPC,indLastNumbMPC+1))>=0) // YES ! 
		{
		  memcpy(myAst[ind].neaSourceRegion, ast.neaSourceRegion, 10*sizeof(ast.neaSourceRegion[0]));
		  //		  myAst[ind].refs[3]+=MASIERO11;
		  cout << ">>match 1  " << ast.pack << "  "<< ind << endl;
		}
	      else
		{ // the designation could not be matched.
		  // try to see if the asteroid has been numbered LATER
		  ast.findId(ids,nIds);
		  if (ast.Number>0)
		    {
		      memcpy(myAst[ast.Number-1].neaSourceRegion, ast.neaSourceRegion, 10*sizeof(ast.neaSourceRegion[0]));
		      //		      myAst[ast.Number-1].refs[3]+=MASIERO11;
		      cout << ">>match 2  " << ast.pack << endl;
		    }
		  else
		    {
		      if ((ind = findUnnumberedAst(myAst,ast,indLastMPC,indLastNumbMPC+1))>=0) // YES ! 
			{
			  memcpy(myAst[ind].neaSourceRegion, ast.neaSourceRegion, 10*sizeof(ast.neaSourceRegion[0]));
			  //			  myAst[ind].refs[3]+=MASIERO11;
			  cout << ">>match 23  " << ast.pack << endl;
			}
		      else
			{
			  cout << ">>match FAILED  " << ast.pack << "--" << line << endl;
			}
		    }
		}
	    }
	}
    myfile.close();
  } else cout << "file nea.source.aei.dat.idl not found " <<endl;
}

template<class T>
void readInMarcoHCMfamily(string filename, long famNumber, T myAst[], long indLastMPC, long indLastNumbMPC, string* ids, long nIds)
{
  ifstream myfile (filename);
  long ind;
  string line;
  AsteroidWithPole ast;
  
  if (myfile.is_open())
  {    
    while ( getline (myfile,line))
      if ((line.length()>1))
	{
	  sscanf(&line.at(0),"%*ld %ld",&ind);
	  myAst[ind-1].familyAdditional=famNumber;
	}
    myfile.close();
  } else cout << "readInMarcoHCMfamily file not found " <<endl;
}

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////// For the moment it stores the thermal inertia in the variable DENSITY and MASS//////////////////////////////////
//////// This is bullshit techinque at the moment.
///////  i read the file and i store the sum of the thermal inertia in the density variable and the sum of the weigts in mass
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template<class T>
void readInThermalInertia(string filename,T myAst[], long indLastMPC, long indLastNumbMPC, string* ids, long nIds)
{
  ifstream myfile (filename);
  long ind;
  string line;
  AsteroidWithPole ast;
  double Gamma, GammaUnc;
  
  if (myfile.is_open())
  {    
    while (getline (myfile,line))
      if ((line.length()>1))
	{
	  sscanf(&line.at(0),"%ld %*s %*d %lf %lf",&ind,&Gamma,&GammaUnc);
	  myAst[ind-1].density+=Gamma*(1./GammaUnc/GammaUnc); // sum the value with the weigth
	  myAst[ind-1].mass+=(1./GammaUnc/GammaUnc); // sum the weigts
	}
    myfile.close();
  } else cout << "readInThermalInertia file not found " <<endl;
}

template<class T>
void readInThermalInertiaHistorical(string filename,T myAst[], long indLastMPC, long indLastNumbMPC, string* ids, long nIds)
{
  ifstream myfile (filename);
  long ind;
  string line;
  AsteroidWithPole ast;
  double Gamma, GammaUnc;
  
  if (myfile.is_open())
  {    
    while (getline (myfile,line))
      if ((line.length()>1))
	{
	  sscanf(&line.at(0),"%ld %*s %*f %*f %lf %lf",&ind,&Gamma,&GammaUnc);
	  if (ind<=0) continue;
	  myAst[ind-1].density+=Gamma*(1./GammaUnc/GammaUnc); // sum the value with the weigth
	  myAst[ind-1].mass+=(1./GammaUnc/GammaUnc); // sum the weigts
	}
    myfile.close();
  } else cout << "readInThermalInertiaHistorical file not found " <<endl;
}

main(int argc, char* argv[])
{

  char filename[]="MPCORB.DAT";
  long N=800000;
  long indLastMPC=0, indLastNumbMPC=0;

  long nIds;
  string* ids=getIds(&nIds);  

  // count the asteroid in the MPCORB.DAT
  countMPCorb(&N,filename);
  cout << "there are " << N << " asteroids" << endl;  
  AsteroidWithPole* myAst = new AsteroidWithPole[N];
  //  AsteroidPelem* myAst = new AsteroidPelem[N];
  //Asteroid* myAst = new Asteroid[N];
  
  readMPCorb(myAst, &indLastMPC, &indLastNumbMPC, filename);
  //  readInPropSYN(myAst);
  cout << indLastMPC << endl;
  cout << indLastNumbMPC << endl;
  cout << indLastMPC-indLastNumbMPC << endl;
  
  cout << " Sorting " << endl; cout.flush();
  // see:
  // http://en.cppreference.com/w/cpp/algorithm/qsort
  // compile with 
  std::qsort(myAst, N, sizeof *myAst, [](const void* a, const void* b)
	     {
	       Asteroid arg1 = *static_cast<const Asteroid*>(a);
	       Asteroid arg2 = *static_cast<const Asteroid*>(b);
	       
	       if(arg1 < arg2) return -1;
	       if(arg1 > arg2) return 1;
	       return 0;
 
	       //  return (arg1 > arg2) - (arg1 < arg2); // possible shortcut
	       //  return arg1 - arg2; // erroneous shortcut (fails if INT_MIN is present)
	     });
  cout << " Done" << endl;

  cout << "Saving data "; cout.flush();
  ast2bin(myAst,N);
  cout << "Done" << endl;
  
  // proper elements
  cout << " Importing Numbered PROPER Elem " << endl; cout.flush();
  readInAllNumANA(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl;

  cout << " Importing UnNumbered PROPER Elem " << endl; cout.flush();
  readInUfitANA(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl;
  
  cout << " Importing Syn pElem " << endl; cout.flush();
  readInAllSYN(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl;

  cout << "Saving data "; cout.flush();
  ast2bin(myAst,N);
  cout << "Done" << endl;
    
  // Tagging asteroids observed by WISE i.e. those that have C51 Observations 
  cout << " Tagging asteroids with C51 reported to MPC " << endl; cout.flush();
  tagAsteroidsWithC51(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();

// BETTER TO SAVE THIS STEP
  cout << "Saving data "; cout.flush();  ast2bin(myAst,N); cout << "Done" << endl;
  
  /// RESTORING DATA !!!!! 
//  cout << "Restoring data "; cout.flush(); AsteroidWithPole* myAst = bin2ast(&N,&indLastMPC,&indLastNumbMPC);  cout << "Done" << endl;

  // clean Diameter and albedo DATA.
  for (long i=0; i<N-1; i++){
    myAst[i].D=-1;
    myAst[i].sigmaD=-1;
    myAst[i].pV=-1;
    myAst[i].sigmapV=-1;
    myAst[i].eta=-1;
    myAst[i].sigmaEta=-1;
    myAst[i].pIR=-1;
    myAst[i].sigmapIR=-1;
    myAst[i].nW1=0;
    myAst[i].nW2=0;
    myAst[i].nW3=0;
    myAst[i].nW4=0;
    myAst[i].refs[3]=0;
  }
  
  // diameters and albedos
  cout << " Importing Nugent 16 " << endl; cout.flush();
  readInNugent16MBA(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  readInNugent16NEA(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();

  cout << " Importing Nugent 15 " << endl; cout.flush();
  readInNugent15(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();
  
  cout << " Importing Ali-Lagoa MarsCrossers " << endl; cout.flush();
  readInVictorMarsCrosser2016(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();

  cout << " Importing Masiero 12 " << endl; cout.flush();
  readInMasiero12(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();
  
  cout << " Importing SIMPS " << endl; cout.flush();
  readInSIMPS(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl;

  cout << " Importing AKARI " << endl; cout.flush();
  readInAkari(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl;
  
  cout << " Importing Masiero 11 " << endl; cout.flush();
  readInMasiero11(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();

  cout << " Importing Masiero 14 " << endl; cout.flush();
  readInMasiero14(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();

  cout << " Importing Ali-Lagoa and Delbo 16 " << endl; cout.flush();
  readInVictor16(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();
  
  cout << "Reading  Earn Table "; cout.flush();
  readInEarnDLR(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << "Done" << endl;
    
  //  cout << " Importing Ali-Lagoa MarsCrossers " << endl; cout.flush();
  //  readInVictorMarsX(myAst,indLastMPC,indLastNumbMPC,ids,nIds); cout << " Done" << endl; cout.flush();

  // BETTER TO SAVE THIS STEP
  cout << "Saving data "; cout.flush();  ast2bin(myAst,N); cout << "Done" << endl;

  // taxonomy
  cout << " Importing Taxonomy 2010" << endl; cout.flush();
  readInTaxonomy10(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();

  cout << " Importing Carvano 2010 taxonomy" << endl; cout.flush();
  readInCarvano10(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();

  cout << " Importing DeMeo & Carry Taxonomy" << endl; cout.flush();
  readInDeMeoCarry13(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();

  cout << " Importing PinillaAlonso+ 2016 Taxonomy" << endl; cout.flush();
  //readInPinillaAlonso16(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();

  cout << " Importing DeLeon+ 2016 Taxonomy" << endl; cout.flush();
  readInDeLeon16(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();
 
  // FIX TAX STRING
  cout << "Fixing Taxonomy data "; cout.flush();
  for (long i=0; i<N; i++)
    if(myAst[i].hasTax)
      myAst[i].fixTaxData();
  cout << "Done" << endl;

  // BETTER TO SAVE THIS STEP
  cout << "Saving data "; cout.flush();  ast2bin(myAst,N); cout << "Done" << endl;
  
  cout << " Importing MOC4  " << endl; cout.flush();
  readInSDSSmoc4(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();
  
  // ROTATION PERIOD
  cout << " Importing LC_SUM_PUB" << endl; cout.flush();
  readInRotationPeriods(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();
  
  // Reset family membership
  cout << " Resetting Devid N family Membership" << endl; cout.flush();
  for (long i=0; i<indLastMPC; i++)
    myAst[i].familyDavidN=0;
  cout << " Done" << endl; cout.flush();
  
  // import David Nesvorny families
  cout << " Importing David N Familes" << endl; cout.flush();
  readInFamiliesDavidN(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();

  /*
    /// RESTORING DATA !!!!! 
  cout << "Restoring data "; cout.flush(); AsteroidWithPole* myAst = bin2ast(&N,&indLastMPC,&indLastNumbMPC);  cout << "Done" << endl;
  */
  
  //  readInMarcoHCMfamily("fam100.list", 8, myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Importing Marco Famiies as David N familes " << endl; cout.flush();
  readInMarcoHCMfamily("fam_84_120.list", 84, myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  readInMarcoHCMfamily("fam_142_140.list", 142, myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  readInMarcoHCMfamily("fam_163_140.list", 163, myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  readInMarcoHCMfamily("fam_302_120.list", 302, myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  readInMarcoHCMfamily("fam_313_110.list", 302, myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  readInMarcoHCMfamily("fam_329_200.list", 329, myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  readInMarcoHCMfamily("fam_495_140.list", 495, myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  readInMarcoHCMfamily("fam_623_230.list", 623, myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  readInMarcoHCMfamily("fam_752_190.list", 752, myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();  
  
  /*
   // read the Olivine abundnace i.e. [ol/(ol+low-Ca px)] %
  cout << "Reading  Olivine abundnace "; cout.flush();
  readInFo(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << "Done" << endl;
  */
  
  /*
  cout << " Importing NEA SourceRegions (Bottke 2005 Model)" << endl; cout.flush();
  readInNEASourceRegion(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();
  */

  /*
  cout << " Resetting NEAsource Regions" << endl; cout.flush();
  for (long i=0; i<indLastMPC; i++)
    myAst[i].resetNEAsSourceRegions();
  cout << " Done" << endl; cout.flush();
  
  cout << " Importing NEA SourceRegions (Granvik+ 2015 model; calc by Morby) " << endl; cout.flush();
  readInNEASourceRegionGranvik(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();
  // the regions in sequence are:
  // nu6; 5/2; 2/1; Hung; 3/1; PHO; JFC
  */
  
  /*
  cout << " readInCibulkova2016 " << endl; cout.flush();
  readInCibulkova2016(myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();
  */

  /*
  cout << " Importing Josef Hanus 2017 thermal inertias " << endl; cout.flush();
  readInThermalInertia("shapes_errorsCor3", myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();  
  cout << " Importing Historical thermal inertias " << endl; cout.flush();
  readInThermalInertiaHistorical("literature6.dat", myAst,indLastMPC,indLastNumbMPC,ids,nIds);
  cout << " Done" << endl; cout.flush();  

  */
  
  cout << "Saving data "; cout.flush(); ast2bin(myAst,N);  cout << "Done" << endl;
 
  for (long i=0; i<10; i++) myAst[i].print();

  return 0;
}
